function [cint] = chebint(a, b, c, n)
% given a, b, and c[1..n] as output from chebft or chebftd, and given n,
% the desired degree of approximation (length of c to be used),
% this routine computes cint, the Chebyshev coefficients of the
% integral of the function whose coeffs are in c. The constant of
% integration is set fo that the integral vanishes at a.
%
sum = 0;
fac = 1.;
con = 0.25*(b-a); % factor that normalizes the interval
for j=2:n-1
   cint(j)=con*(c(j-1)-c(j+1))/(j-1);
   sum = sum + fac * cint(j);
   fac = - fac;
end
cint(n) = con*c(n-1)/n;
sum = sum + fac*cint(n);
cint(1) = 2.0*sum; % set constant of integration.
return;
